mask.dir = '/home1/99/jc152199/ESA/MASK/'

# Define Directory of ASCII files

BC.dir = "/home1/99/jc152199/ESA/images/"

# Define the out directory
 
out.dir = "/home1/99/jc152199/ESA/images/"

# Read in the template ASCII (projection is UTM WGS1984 Projected (i.e. easting/northing)

base.asc = read.asc(paste(mask.dir,'mask.asc',sep=""))

# Convert all values to zero 

base.asc = base.asc * 0

# Convert base.asc into a dataframe with 1 row for each row/column position

base.pos = as.data.frame(which(is.finite(base.asc), arr.ind = T)) 

# Add east and north to this data frame of row/column positions

base.pos$east = getXYcoords(base.asc)$x[base.pos$row]
base.pos$north = getXYcoords(base.asc)$y[base.pos$col]

# Read in an ASCII file to be resized

t.asc = read.asc(paste(BC.dir,layers[1],sep=""))

# Create a two column data frame of row/column positions

t.pos =  as.data.frame(which(is.finite(t.asc), arr.ind = T))

# Calculate east/north for all row/column positions
    
t.pos$east = getXYcoords(t.asc)$x[t.pos$row]
t.pos$north = getXYcoords(t.asc)$y[t.pos$col]

# List files to resize & project

layers = list.files(BC.dir, pattern = '.asc')

  for (layer in layers) 
  
    {
    
    t.asc = read.asc(paste(BC.dir,layer,sep=""))
    
    t.data = extract.data(cbind(base.pos$east,base.pos$north),t.asc)
    
    cat(head(layer),'...\n')

    base.asc[cbind(base.pos$row,base.pos$col)] = t.data
            
    write.asc(base.asc, paste(out.dir,layer,sep=""))
    
    cat(head(layer),'...\n')
    
    } 

################################################################################
# This script was created to recreate microclimate based on a linear model
# This script will read in some static ASCII data (coastdist, EVI)
# And loop through Tmax and Solar data from AWAP

################################################################################
#load necessary libraries

#list the libraries needed
necessary=c("adehabitat","SDMTools","rgdal","sp")

#check if library is installed
installed = necessary %in% installed.packages()

#if library is not installed, install it
if (length(necessary[!installed]) >=1) install.packages(necessary[!installed], dep = T)

#load the libraries
for (lib in necessary) library(lib,character.only=T)

################################################################################

# Define Directory of ASCII files

ASCII.dir = "/home1/99/jc152199/ESA/"

# Define the out directory
 
out.dir = "/home1/99/jc152199/ESA/images/"

# Create list of layers

layers = list.files(ASCII.dir, pattern='.asc')

# Read in ASCII files

AWAPdaily.asc = read.asc(paste(ASCII.dir,layers[1],sep=""))
cd.asc = read.asc(paste(ASCII.dir,layers[2],sep=""))
EVI.asc = read.asc(paste(ASCII.dir,layers[3],sep=""))
solar.asc = read.asc(paste(ASCII.dir,layers[4],sep=""))
total.asc = read.asc(paste(ASCII.dir,layers[6],sep=""))
output.asc = read.asc(paste(ASCII.dir,'output.asc',sep=""))
totalsmall.asc = read.asc(paste(ASCII.dir,'Resize/','above28cuonly.asc',sep=""))

# Define an ASCII to represent Under Log Temp as AWAPdaily.asc - 6

ULdaily.asc = output.asc -6    

                              
png(paste(out.dir,'daysabove28cuonly.png',sep=""),width=600,height=1000,bg=rgb(.5,.5,.5,0,names='clear')) #Command bg=rgb(.5,.5,.5,0,names='clear') produces a transparent colour
      
      par(pty="m",oma=c(5,0,5,0),mar=c(0,0,0,0),cex=1,xpd=T)
      
      image(totalsmall.asc,axes=F,ann=F,frame.plot=F,oldstyle =T,col=c(heat.colors(9)[9:1]))
      
      cat(layers[i],'...\n')
      
      dev.off()
      
png(paste(out.dir,'ULDailyTmax.png',sep=""),width=600,height=1000,bg=rgb(.5,.5,.5,0,names='clear')) #Command bg=rgb(.5,.5,.5,0,names='clear') produces a transparent colour
      
      par(pty="m",oma=c(5,0,5,0),mar=c(0,0,0,0),cex=1,xpd=T)
      
      image(ULdaily.asc,axes=F,ann=F,frame.plot=F,breaks=c(17:36),oldstyle =T,col=c(heat.colors(19)[19:1]))
      
      cat(layers[i],'...\n')
      
      dev.off()
      
png(paste(out.dir,'AirDailyTmax.png',sep=""),width=600,height=1000,bg=rgb(.5,.5,.5,0,names='clear')) #Command bg=rgb(.5,.5,.5,0,names='clear') produces a transparent colour
      
      par(pty="m",oma=c(5,0,5,0),mar=c(0,0,0,0),cex=1,xpd=T)
      
      image(output.asc,axes=F,ann=F,frame.plot=F,breaks=c(17:36),oldstyle =T,col=c(heat.colors(19)[19:1]))
      
      cat(layers[i],'...\n')
      
      dev.off()
      
png(paste(out.dir,'cd.png',sep=""),width=600,height=1000,bg=rgb(.5,.5,.5,0,names='clear')) #Command bg=rgb(.5,.5,.5,0,names='clear') produces a transparent colour
      
      par(pty="m",oma=c(5,0,5,0),mar=c(0,0,0,0),cex=1,xpd=T)
                                                     
      image(cd.asc,axes=F,ann=F,frame.plot=F,oldstyle =T,col=c(grey(c(0,.1,.2,.3,.4,.5,.6,.7,.8,.9,1)[1:11])))
      
      dev.off()     

png(paste(out.dir,'solar.png',sep=""),width=600,height=1000,bg=rgb(.5,.5,.5,0,names='clear')) #Command bg=rgb(.5,.5,.5,0,names='clear') produces a transparent colour
      
      par(pty="m",oma=c(5,0,5,0),mar=c(0,0,0,0),cex=1,xpd=T)
      
      image(solar.asc,axes=F,ann=F,frame.plot=F,oldstyle =T,col=c(grey(c(0,.1,.2,.3,.4,.5,.6,.7,.8,.9,1)[1:11])))
      
      dev.off()
      
png(paste(out.dir,'total.png',sep=""),width=600,height=1000,bg=rgb(.5,.5,.5,0,names='clear')) #Command bg=rgb(.5,.5,.5,0,names='clear') produces a transparent colour
      
      par(pty="m",oma=c(5,0,5,0),mar=c(0,0,0,0),cex=1,xpd=T)
      
      image(total.asc,axes=F,ann=F,frame.plot=F,oldstyle =T,col=c(heat.colors(16)[16:1]))
      
      dev.off()
      
png(paste(out.dir,'EVI.png',sep=""),width=600,height=1000,bg=rgb(.5,.5,.5,0,names='clear')) #Command bg=rgb(.5,.5,.5,0,names='clear') produces a transparent colour
      
      par(pty="m",oma=c(5,0,5,0),mar=c(0,0,0,0),cex=1,xpd=T)
      
      image(EVI.asc,axes=F,ann=F,frame.plot=F,oldstyle =T,col=c(cols))
      
      dev.off()
      
# Create color pallete for EVI map      
      
greens = colors()[t.vec]
t.vec = c(1,18,20,83,88,85,81,47)
d=col2rgb(greens[2:7])

z=rgb(d[1,1], d[2,1], d[3,1], maxColorValue=255)
y=rgb(d[1,2], d[2,2], d[3,2], maxColorValue=255)
x=rgb(d[1,3], d[2,3], d[3,3], maxColorValue=255)
w=rgb(d[1,4], d[2,4], d[3,4], maxColorValue=255)
v=rgb(d[1,5], d[2,5], d[3,5], maxColorValue=255)
u=rgb(d[1,6], d[2,6], d[3,6], maxColorValue=255)
uu=rgb(d[1,7], d[2,7], d[3,7], maxColorValue=255)

cols = c(z,y,x,w,v,u,uu)



d= col2rgb(c("yellow","red"))

